function gs_strashow(GSS,visScale)
%gs_strashow - strain field visualization
%
%       gs_strashow(GSS,visScale)
%
% This function displays the strain field vectors included in GSS struct
% variable on the basis of the visualization scale visScale.
% If visScale is undefined or empty, it is 
%
%   visScale = min(dx,dy)/max(mean(abs(EMAX(:))),mean(abs(EMIN(:)))),
%
% where all parameters are taken from GSS. 
%
% These plots are shown:
%   1. significance areas;
%   2. strain field (in all grid nodes, high significance grid nodes or
%      high and mean significance grid nodes, on the basis of user choice.
%      This choice also applies to plots 3-6 and 8);
%   3. change in area contour map;
%   4. prevailing eigenvaue contour map;
%   5. normalized shear contour map;
%   6. second invariant of the strain;
%   7. error on maximum and minimum strain contour maps;
%   8. strain field in high/high and mean significance area overlapped on
%      prevailing eigenvalue contour map
%   9. velocity/translation field.
%
% See also GS_ITERA, GS_STRAIN.

% G. Teza & A. Pesci, 2004, 2005, 2008, 2011, 2022.

disp('Visualization in progress...');

Options = GridStrainOptions;

LW = 1.5;

% FontSize single visualizations:
FontSizeLabel = 16;     % FontSize labels
FontSizeTitle = 16;     % FontSize titles
FontSizeAxis = 14;      % FontSize axis

e = GSS.East;
n = GSS.North;
ve = GSS.EastVelocity;
vn = GSS.NorthVelocity;
IMS = GSS.RefImageStruct;
npe = GSS.IncludedPoints;
d0 = GSS.ScaleFactor;
EMAX = GSS.MaxStrainMatrix;
EMIN = GSS.MinStrainMatrix;
PHI = GSS.AzimuthStrainMatrix;
eEMAX = GSS.MaxStrainMatrixError;
eEMIN = GSS.MinStrainMatrixError;
FLC = GSS.ComputationFlag;
FLS = GSS.SignificanceFlag;
U0 = GSS.UTranslation;
V0 = GSS.VTranslation;
eU0 = GSS.UTranslationError;
eV0 = GSS.VTranslationError;
OMEG = GSS.Rotation;
eOMEG = GSS.RotationError;
dx = GSS.XGridStep;
dy = GSS.YGridStep;
ns = GSS.StationID;
stn = GSS.StationName;

ei = e; 
ni = n;
vei = ve;
vni = vn;
ei(~npe) = NaN; 
ni(~npe) = NaN;
vei(~npe) = NaN; 
vni(~npe) = NaN;
en = e; 
nn = n;
ven = ve;
vnn = vn;
en(npe) = NaN; 
nn(npe) = NaN;
ven(npe) = NaN; 
vnn(npe) = NaN;

if ~isempty(ns)
    nsi = ns;
    nsi(~npe) = NaN;
    nsn = ns;
    nsn(npe) = NaN;
    stni = [];
    stnn = [];
else
    nsi = [];
    nsn = [];
    stni = stn;
    stni(~npe) = {[]};
    stnn = stn;
    stnn(npe) = {[]};
end

[X,Y] = gs_gridExtra(GSS);
szx = size(X);

EIG1B = nan(szx(1),szx(2),2);   % In order to draw the first eigenvector
EIG1R = nan(szx(1),szx(2),2);   % B (blue): dilatation; R (red): compression   
EIG2B = nan(szx(1),szx(2),2);
EIG2R = nan(szx(1),szx(2),2); 
for m = 1:szx(1)
    for mm = 1:szx(2)
        if FLC(m,mm) == 1
            if PHI(m,mm) >= 0
                phi = PHI(m,mm);
            else
                phi = pi+PHI(m,mm);
            end
            aema = abs(EMAX(m,mm));
            aemi = abs(EMIN(m,mm));
            cphi = cos(phi);
            sphi = sin(phi);
            if phi <= pi/2   % Angle in the range 0-90�
                if EMAX(m,mm) >= 0
                    EIG1B(m,mm,1) = aema*sphi;
                    EIG1B(m,mm,2) = aema*cphi;
                else
                    EIG1R(m,mm,1) = aema*sphi;
                    EIG1R(m,mm,2) = aema*cphi;
                end
                if EMIN(m,mm) >= 0
                    EIG2B(m,mm,1) = aemi*cphi;
                    EIG2B(m,mm,2) = -aemi*sphi;
                else
                    EIG2R(m,mm,1) = aemi*cphi;
                    EIG2R(m,mm,2) = -aemi*sphi;
                end
            else
                if EMAX(m,mm) >= 0
                    EIG1B(m,mm,1) = aema*sphi;  % sin(pi-phi) = sin(phi)
                    EIG1B(m,mm,2) = aema*cphi;  % cos(pi-phi) = -cos(phi) 
                else
                    EIG1R(m,mm,1) = aema*sphi;
                    EIG1R(m,mm,2) = aema*cphi;
                end
                if EMIN(m,mm) >= 0
                    EIG2B(m,mm,1) = -aemi*cphi; 
                    EIG2B(m,mm,2) = aemi*sphi;
                else
                    EIG2R(m,mm,1) = -aemi*cphi;
                    EIG2R(m,mm,2) = aemi*sphi;
                end
            end
        end
    end
end

figOption = figure('Name','Visualization option');
set(figOption,'Units','normalized','OuterPosition',[0.1 0.1 0.7 0.5]);
imshow([])
title('VISUALIZATION OPTION','FontSize',24,'Color','k');  
axis equal;
lsOption = {...
    'PLOTS IN ALL GRID NODES',...
    'PLOTS IN HIGH GEOMETRIC SIGNIFICANCE GRID NODES',...
    'PLOTS IN HIGH AND MEAN GEOMETRIC SIGNIFICANCE GRID NODES',...
    'NO PLOTS'};
lbOption = uicontrol(figOption,'Style','listbox',...
    'string',lsOption,...
    'FontSize',14,...
    'Units','normalized','Position', [0.2 0.3 0.60 0.40],...
    'Callback', 'uiresume(gcbf)');
uiwait(figOption);
optFig = lbOption.Value;
delete(lbOption);
close(figOption);

if optFig == 1
    optStr1 = 'all';
    optStr2 = '';
    strAHM = 'ALL GRID NODES';
    IAHM = true(size(EMAX));
elseif optFig == 2
    optStr1 = 'high';
    optStr2 = '';
    strAHM = 'HIGH SIGNIFICANCE GRID NODES';
    IAHM = FLS == 2; 
elseif optFig == 3
    optStr1 = 'high';
    optStr2 = 'medium';
    strAHM = 'HIGH AND MID SIGNIFICANCE GRID NODES';
    IAHM = FLS > 0;
else
    return
end
maemax = max([abs(EMAX(IAHM(:)));abs(EMIN(IAHM(:)))]);
maemax = ApprInt(maemax);

if nargin < 2
    visScale = [];
end
if isempty(visScale)
    meabsma  = max(abs(EMAX(IAHM(:))));
    meabsmi  = max(abs(EMIN(IAHM(:))));
    visScale = min(dx,dy)/max(meabsma,meabsmi);
end

XAHM = X;
XAHM(~IAHM) = NaN;
YAHM = Y;
YAHM(~IAHM) = NaN;

%% Velocity vectors (only if at least a point is excluded from computations)

if sum(npe) < numel(npe)
    ALA0 = gs_velPlot(GSS);
    hold on;
    gs_ReadBGImage(IMS,ALA0); 
    hold off
end

%% SIGNIFICANCE MAP

whitmap = [1 1 1];               % RGB white
graymap = [.5 .5 .5];            % gray
aquamap = [127/255 1 212/255];   % aquamarine 
mapSign = [whitmap; graymap; aquamap];
figSign = figure; 
set(figSign,'Units','normalized','OuterPosition',[0 0 1 1]);  
hold on;
% contourf(X,Y,FLS,3); 
surf(X,Y,FLS);
colormap(mapSign); 
colorbar;
plot(X,Y,'g',X',Y','g');
if Options.VisualizePoints
    plot(ei,ni,'xk',en,nn,'xm'); 
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,FontSizeAxis);
    end
end
title(sprintf('%s\n%s%s%s',...
    'SIGNIFICANCE OF CALCULATIONS ON GRID NODES',...
    'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',FontSizeTitle);
xlabel('EAST (m)','FontSize',FontSizeLabel); 
ylabel('NORTH (m)','FontSize',FontSizeLabel); 
axis equal; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; 
set(gca,'FontSize',12);

%% STRAIN/STRAIN RATE 

figs = figure; 
set(figs,'Units','normalized','OuterPosition',[0 0 1 1]); 
hold on;

if Options.VisualizePoints
    ppl = plot(ei,ni,'xk',en,nn,'xm');
    legend([ppl(1) ppl(2)],'Point included','Point not included',...
        'Location','NorthEastOutside');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,FontSizeAxis);
    end
end
    
% plot of strain field by means of ancillary function (all or high):
PlotStrainField(X,Y,EIG1B,EIG1R,EIG2B,EIG2R,FLS,visScale,Options,optStr1);
% medium (not performed if optStr2 is empty)
PlotStrainField(X,Y,EIG1B,EIG1R,EIG2B,EIG2R,FLS,visScale,Options,optStr2);

%--- units representation:

quiver(X(szx(1),szx(2))+2*dx,Y(1,szx(2)),visScale*maemax,0,0,'b',...
    'HandleVisibility','off');
text(X(szx(1),szx(2))+dx,Y(1,szx(2))+dy/2,...
    sprintf('STRAIN SCALE: %.3g 1/y',maemax));
xlabel('EAST (m)','FontSize',FontSizeLabel); 
ylabel('NORTH (m)','FontSize',FontSizeLabel); 
title(sprintf('%s\n%s%s%s%s',...
    'STRAIN FIELD (BLUE: EXTENSION, RED: COMPRESSION)',...
    strAHM,' - SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',FontSizeTitle);
axis equal; 

ALA = axis;
gs_ReadBGImage(IMS,ALA);
set(gca,'FontSize',FontSizeAxis);
hold off; 

%% CHANGE IN AREA CONTOUR MAP   

figca = figure; 
set(figca,'Units','normalized','OuterPosition',[0 0 1 1]); 
hold on;
TrE = EMAX+EMIN;
TrE(~IAHM) = NaN;
contourMod(X,Y,TrE,Options);
colorbar; 
colormap(flipud(jet)); 

if Options.VisualizePoints
    plot(ei,ni,'xk',en,nn,'xm');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,FontSizeAxis);
    end
end
if (optFig == 1) || ((optFig > 1) && Options.VisualizeGridHSA)
    plot(XAHM,YAHM,'g',XAHM',YAHM','g');
end

xlabel('EAST (m)','FontSize',FontSizeLabel);
ylabel('NORTH (m)','FontSize',FontSizeLabel);
axis image; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; set(gca,'FontSize',FontSizeAxis);
title(sprintf('%s - %s\n%s%s%s',...
    'CHANGE IN AREA CONTOUR MAP',strAHM,...
    'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',FontSizeTitle);

%% PREVAILING EIGENVALUE FOR EACH GRID NODE (WITH ITS SIGN) 

figpe = figure; 
set(figpe,'Units','normalized','OuterPosition',[0 0 1 1]); 
hold on;
ME = nan(size(X));
for m = 1:szx(1)
    for mm = 1:szx(2)
        if FLC(m,mm) == 1
            if EMAX(m,mm) >= 0
                if EMAX(m,mm) >= abs(EMIN(m,mm))
                    ME(m,mm) = EMAX(m,mm);
                else
                    ME(m,mm) = EMIN(m,mm);
                end
            else
                ME(m,mm) = EMIN(m,mm);
            end
        end
    end
end
ME(~IAHM) = NaN;
contourMod(X,Y,ME,Options);
colorbar; 
colormap(flipud(jet)); 

if (optFig == 1) || ((optFig > 1) && Options.VisualizeGridHSA)
    plot(XAHM,YAHM,'g',XAHM',YAHM','g');
end
if Options.VisualizePoints
    plot(ei,ni,'xk',en,nn,'xm');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,FontSizeAxis);
    end
end
xlabel('EAST (m)','FontSize',FontSizeLabel); 
ylabel('NORTH (m)','FontSize',FontSizeLabel);
axis image; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; set(gca,'FontSize',FontSizeAxis);
title(sprintf('%s - %s\n%s%s%s',...
    'PREVAILING EIGENVALUE',strAHM,...
    'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',FontSizeTitle);

%% NORMALIZED SHEAR CONTOUR MAP

figins = figure; 
set(figins,'Units','normalized','OuterPosition',[0 0 1 1]); 
hold on
ShE = EMAX-EMIN;
Inshe = abs(TrE) < 10*eps;
Ishe  = ~Inshe;
ShE(Ishe)  = ShE(Ishe)./TrE(Ishe);
ShE(Inshe) = NaN;
ShE(~IAHM) = NaN;
contourMod(X,Y,ShE,Options);
colorbar; 
colormap(flipud(jet));

if (optFig == 1) || ((optFig > 1) && Options.VisualizeGridHSA)
    plot(XAHM,YAHM,'g',XAHM',YAHM','g');
end
if Options.VisualizePoints
    plot(ei,ni,'xk',en,nn,'xm');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,FontSizeAxis);
    end
end
xlabel('EAST (m)','FontSize',FontSizeLabel); 
ylabel('NORTH (m)','FontSize',FontSizeLabel);
axis image; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; set(gca,'FontSize',FontSizeAxis);
title(sprintf('%s - %s\n%s%s%s',...
    'NORMALIZED SHEAR CONTOUR MAP',strAHM,...
    'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',FontSizeTitle);


%% SECOND INVARIANT OF THE STRAIN

figca = figure; 
set(figca,'Units','normalized','OuterPosition',[0 0 1 1]); 
hold on;
SIS = sqrt(EMAX.^2+EMIN.^2);
SIS(~IAHM) = NaN;
contourMod(X,Y,SIS,Options);
colorbar; 
colormap(flipud(jet));  

if Options.VisualizePoints
    plot(ei,ni,'xk',en,nn,'xm');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,FontSizeAxis);
    end
end
if (optFig == 1) || ((optFig > 1) && Options.VisualizeGridHSA)
    plot(XAHM,YAHM,'g',XAHM',YAHM','g');
end

xlabel('EAST (m)','FontSize',FontSizeLabel);
ylabel('NORTH (m)','FontSize',FontSizeLabel);
axis image; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; set(gca,'FontSize',FontSizeAxis);
title(sprintf('%s - %s\n%s%s%s',...
    'SECOND INVARIANT OF THE STRAIN',strAHM,...
    'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',FontSizeTitle);

%% RELATIVE ERROR CONTOUR MAPS 

% maximum strain:
figrec = figure; 
set(figrec,'Units','normalized','OuterPosition',[0 0 1 1]);  
subplot(1,2,1);
hold on;
Inan   = isnan(EMAX);
Ienan  = isnan(eEMAX);
Inode  = abs(EMAX) < 10*eps;
Inosi  = FLS == 0;
Inobo  = logical(Inan+Ienan+Inode+Inosi);
Ishe   = ~Inobo;
EREMA  = nan(szx);
EREMA(Ishe) = log10(eEMAX(Ishe)./abs(EMAX(Ishe)));
contourMod(X,Y,real(EREMA),Options);
colorbar; 
colormap jet;

if (optFig == 1) || ((optFig > 1) && Options.VisualizeGridHSA)
    plot(XAHM,YAHM,'g',XAHM',YAHM','g');
end
if Options.VisualizePoints
    plot(ei,ni,'xk',en,nn,'xm');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,...
            round(0.85*FontSizeAxis));
    end
end

xlabel('EAST (m)','FontSize',round(0.85*FontSizeLabel));
ylabel('NORTH (m)','FontSize',round(0.85*FontSizeLabel));
axis image; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; 
set(gca,'FontSize',round(0.85*FontSizeAxis));
title(sprintf('%s\n%s%s%s',...
    'RELATIVE ERROR (log_{10}(\it\sigma_e\rm_{max}/\ite\rm_{max}))',...
    'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',round(0.85*FontSizeTitle));

% minimum strain:
subplot(1,2,2);
hold on;
Inan   = isnan(EMIN);
Ienan  = isnan(eEMIN);
Inode  = abs(EMIN) < 10*eps;
Inosi  = FLS == 0;
Inobo  = logical(Inan+Ienan+Inode+Inosi);
Ishe   = ~Inobo;
EREMI  = nan(szx);
EREMI(Ishe) = log10(eEMIN(Ishe)./abs(EMIN(Ishe)));
contourMod(X,Y,real(EREMI),Options);
colorbar; 
colormap jet;
if (optFig == 1) || ((optFig > 1) && Options.VisualizeGridHSA)
    plot(XAHM,YAHM,'g',XAHM',YAHM','g');
end
if Options.VisualizePoints
    plot(ei,ni,'xk',en,nn,'xm');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,...
            round(0.85*FontSizeAxis));
    end
end
xlabel('EAST (m)','FontSize',round(0.85*FontSizeLabel)); 
ylabel('NORTH (m)','FontSize',round(0.85*FontSizeLabel));
axis image; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; 
set(gca,'FontSize',round(0.85*FontSizeAxis));
title(sprintf('%s\n%s%s%s',...
    'RELATIVE ERROR (log_{10}(\it\sigma_e\rm_{min}/\ite\rm_{min}))',...
    'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',round(0.85*FontSizeTitle));

%% TRANSLATIONS/VELOCITIES

figvel = figure; 
set(figvel,'Units','normalized','OuterPosition',[0 0 1 1]); 
hold on;
if (optFig == 1) || ((optFig > 1) && Options.VisualizeGridHSA)
    plot(XAHM,YAHM,'g',XAHM',YAHM','g');
end
U0AHM = U0;
U0AHM(~IAHM) = NaN;
V0AHM = V0;
V0AHM(~IAHM) = NaN;
quiver(X,Y,U0AHM,V0AHM,'HandleVisibility','off');
if Options.VisualizePoints
    ppv = plot(ei,ni,'xk',en,nn,'xm');
    legend([ppv(1) ppv(2)],'Point included','Point not included',...
        'Location','NorthEastOutside');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,FontSizeAxis);
    end
end

%--- units representation:
maemaxu = max(abs(U0(:))); 
maemaxv = max(abs(V0(:)));
maemaxuv = max(maemaxu,maemaxv); 
quiver(X(szx(1),szx(2))+2*dx,Y(1,szx(2)),maemaxuv,0,'HandleVisibility','off');
text(X(szx(1),szx(2))+dx,Y(1,szx(2))+dy/2,sprintf('MAX VELOCITY: %.3f mm/y',maemaxuv));
xlabel('EAST (m)','FontSize',FontSizeLabel); 
ylabel('NORTH (m)','FontSize',FontSizeLabel); 
title(sprintf('%s - %s\n%s%s%s',...
    'VELOCITY FIELD',strAHM,...
    'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',FontSizeTitle);
axis equal; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; 
set(gca,'FontSize',FontSizeAxis);

%% PLOT OF STRAIN FIELD AND PREVAILING EIGENVALUE CONTOURING

% High significance only: 

figss = figure; 
set(figss,'Units','normalized','OuterPosition',[0 0 1 1]);

subplot(1,2,1);
hold on;

MEMH = ME;
IH = FLS == 2;
MEMH(~IH) = NaN;
[~,~] = contourf(X,Y,MEMH,10); 
colorbar; 
colormap(flipud(jet)); 
if Options.VisualizePoints
    plot(ei,ni,'xk',en,nn,'xm');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,...
            round(0.85*FontSizeAxis));
    end
end
   
PlotStrainField(X,Y,EIG1B,EIG1R,EIG2B,EIG2R,FLS,visScale,Options,'high');

%--- scale arrow
quiver(X(szx(1),3)+2*dx,Y(3,szx(2)),...
    visScale*maemax,0,0,'b','LineWidth',LW);
text(X(szx(1),3)+dx,Y(3,szx(2))+dy/2,sprintf('STRAIN RATE SCALE: %.3g 1/y',maemax));

title(sprintf('%s - %s\n%s%s%s',...
    'STRAIN FIELD AND PREVAILING EIGENVALUE','HIGH SIGNIFICANCE',...
    'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',round(0.85*FontSizeTitle));
    
xlabel('EAST (m)','FontSize',round(0.85*FontSizeLabel)); 
ylabel('NORTH (m)','FontSize',round(0.85*FontSizeLabel)); 
axis image; 
axis equal; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; 
set(gca,'FontSize',round(0.85*FontSizeAxis));

% High and mean significance: 

subplot(1,2,2);
hold on;

MEMHM = ME;
IHM = FLS > 0;
MEMHM(~IHM) = NaN;
[~,~] = contourf(X,Y,MEMHM,10);
colorbar; 
colormap jet; 
if Options.VisualizePoints
    plot(ei,ni,'xk',en,nn,'xm');
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,en,nn,ven,vnn,nsn,stnn,...
            round(0.85*FontSizeAxis));
    end
end
   
PlotStrainField(X,Y,EIG1B,EIG1R,EIG2B,EIG2R,FLS,visScale,Options,'high');
PlotStrainField(X,Y,EIG1B,EIG1R,EIG2B,EIG2R,FLS,visScale,Options,'medium');    

%--- scale arrow
quiver(X(szx(1),3)+2*dx,Y(3,szx(2)),...
    visScale*maemax,0,0,'b','LineWidth',LW);
text(X(szx(1),3)+dx,Y(3,szx(2))+dy/2,sprintf('STRAIN RATE SCALE: %.3g mm/y',maemax));

title(sprintf('%s - %s\n%s%s%s',...
        'STRAIN FIELD AND PREVAILING EIGENVALUE','HIGH AND MEAN SIGNIFICANCE',...
        'SCALE FACTOR: ',orderapp(d0,6),' m'),'FontSize',round(0.85*FontSizeTitle));
    
xlabel('EAST (m)','FontSize',round(0.85*FontSizeLabel)); 
ylabel('NORTH (m)','FontSize',round(0.85*FontSizeLabel)); 
axis image; 
axis equal; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; 
set(gca,'FontSize',round(0.85*FontSizeAxis));


%% SINGLE POINT OR AREA INTERACTIVE SELECTION
 
xgri = X(1,:);
ygri = Y(:,1);
nosing = 0; 
while nosing == 0
    
    figure(figs);
    if exist('figFin','var')
        figure(figFin); % to have the report (if existing) on foreground
    end
    lbsel = uicontrol(figs,'Style','listbox',...
        'string',...
        {'SHOW THE RESULTS IN A SELECTED GRID POINT',...
            'SHOW MEAN RESULTS IN A SELECTED POLYGONAL AREA',...
            'SKIP'},...
        'FontSize',14,...
        'Units','normalized','Position', [0.02 0.02 0.40 0.20],...
        'Callback', 'uiresume(gcbf)');
    uiwait(figs);
    menselp = lbsel.Value;
    delete(lbsel);

    if menselp == 1
        
        [xpick,ypick] = ginput(1);
        [~,kkmi] = min(abs(xpick-xgri));
        [~,kmi]  = min(abs(ypick-ygri));
        if length(kkmi) > 1
            kkmi = min(kkmi); 
        end
        if length(kmi) > 1
            kmi = min(kmi);   
        end
        if FLC(kmi,kkmi) == 0
            rep_point = sprintf(...
                'GRID-NODE COORDINATES: \n east: %f, \north: %f;\n %s.',...
                X(kmi,kkmi),Y(kmi,kkmi),...
                '=== IN THIS POINT, CALCULATIONS CANNOT BE PERFORMED ===');
        else
            %--------------------------------------------------------------
            rep_point = gs_reppoint(X(kmi,kkmi),Y(kmi,kkmi),...
                EMAX(kmi,kkmi),EMIN(kmi,kkmi),...
                eEMAX(kmi,kkmi),eEMIN(kmi,kkmi),...
                PHI(kmi,kkmi),d0,FLS(kmi,kkmi),...
                U0(kmi,kkmi),V0(kmi,kkmi),eU0(kmi,kkmi),eV0(kmi,kkmi),...
                OMEG(kmi,kkmi),eOMEG(kmi,kkmi));
            %--------------------------------------------------------------
        end
        fprintf('\n----------------------------------------------');
        disp(rep_point);
        fprintf('----------------------------------------------\n');

        figFin = figure('Name','GridStrain');
        set(figFin,'Units','normalized','OuterPosition',[0.1 0.1 0.7 0.7]);
        imshow([]);
        title(rep_point,'FontSize',12,'Color','k');  
    
    elseif menselp == 2
        
        Isel = selePolyPoints(X,Y,figs);
        if sum(Isel(:)) == numel(Isel)
            Isel = Isel(:);
        end
        if (sum(Isel(:)) > 1) && (sum(FLC(Isel(:)) == 1)> 1)
            miX = min(X(Isel(:)));
            maX = max(X(Isel(:)));
            miY = min(Y(Isel(:)));
            maY = max(Y(Isel(:)));
            meEMAX = mean(EMAX(Isel));
            meEMIN = mean(EMIN(Isel));
            sdEMAX = std(EMAX(Isel));
            sdEMIN = std(EMIN(Isel));
            % (a problem on angle must be solved to have meaningful mean)
            mePHI = nanmean(abs(PHI(Isel)))*180/pi;
            npPHI = sum(sign(PHI(Isel)) >= 0);
            nnPHI = sum(sign(PHI(Isel)) < 0);
            if npPHI < nnPHI
                mePHI = -mePHI;
            end
            meTrE = nanmean(TrE(Isel)); %#ok<*NANMEAN> 
            meShE = nanmean(ShE(Isel));
            meSIS = nanmean(SIS(Isel));
            nH = sum(FLS(Isel(:)) == 2);
            nM = sum(FLS(Isel(:)) == 1);
            nL = sum(FLS(Isel(:)) == 0);
            rep_area = gs_reparea(miX,maX,miY,maY,meEMAX,meEMIN,...
                sdEMAX,sdEMIN,mePHI,meTrE,meShE,meSIS,d0,nH,nM,nL);
        else
            rep_area = 'NO AREA DATA CAN BE SHOWN';
        end
        
        fprintf('\n----------------------------------------------');
        disp(rep_area);
        fprintf('----------------------------------------------\n');

        figFin = figure('Name','GridStrain');
        set(figFin,'Units','normalized','OuterPosition',[0.1 0.1 0.7 0.7]);
        imshow([]);
        title(rep_area,'FontSize',12,'Color','k'); 
        
    else
        
        nosing = 1; 
    
    end
    
end

disp('... Visualization completed');

%% ANCILLARY FUNCTIONS

function addIDName(ei,ni,vei,vni,nsi,stni,ex,nx,vex,vnx,nsx,stnx,FSA) 
% Add station ID or station name
mediff = mean([nanstd(diff(ei)) nanstd(diff(ni))]); %#ok<NANSTD>
for kk = 1:numel(ei)
    if ~isempty(nsi) && ~isnan(nsi(kk))
        strikk = num2str(nsi(kk));
        okikk = true;
    elseif ~isempty(stni) && ~isempty(stni{kk})
        strikk = stni{kk};
        okikk = true;
    else
        okikk = false;
    end
    if okikk 
        se = sign(vei(kk)); 
        sn = sign(vni(kk));
        ht = text(ei(kk)-se*mediff/12,ni(kk)-sn*mediff/12,strikk);
        set(ht,'Color','k','FontSize',FSA);
    end
end
for kk = 1:numel(ex)
    if ~isempty(nsx) && ~isnan(nsx(kk))
        strxkk = num2str(nsx(kk));
        okxkk = true;
    elseif ~isempty(stnx) && ~isempty(stnx{kk})
        strxkk = stnx{kk};
        okxkk = true;
    else
        okxkk = false;
    end
    if okxkk 
        se = sign(vex(kk)); sn = sign(vnx(kk));
        ht = text(ex(kk)-se*mediff/12,nx(kk)-sn*mediff/12,strxkk);
        set(ht,'Color','m','FontSize',FSA);
    end
end

function PlotStrainField(X,Y,EIG1B,EIG1R,EIG2B,EIG2R,FLS,visScale,Options,optVis)
% Plot strain vectors
% Initial options management:
LW = 1.5;
switch optVis
    case 'all'
        strLineB = 'b';
        strLineR = 'r';
        optSel = false;
    case 'high'
        strLineB = 'b';
        strLineR = 'r';
        Isel = FLS == 2;
        optSel = true;
    case 'medium'
        strLineB = ':b';
        strLineR = ':r';
        Isel = FLS == 1;
        optSel = true;
    otherwise
        return
end
% Blue and red arrays management:
if optSel
    % Eigenvalue 1, blue:
    EIG1B1 = EIG1B(:,:,1);
    EIG1B2 = EIG1B(:,:,2);
    EIG1B1(~Isel) = NaN;
    EIG1B2(~Isel) = NaN;
    EIG1B(:,:,1) = EIG1B1;
    EIG1B(:,:,2) = EIG1B2;
    % Eigenvalue 2, blue:
    EIG2B1 = EIG2B(:,:,1);
    EIG2B2 = EIG2B(:,:,2);
    EIG2B1(~Isel) = NaN;
    EIG2B2(~Isel) = NaN;
    EIG2B(:,:,1) = EIG2B1;
    EIG2B(:,:,2) = EIG2B2;
    % Eigenvalue 1, red:
    EIG1R1 = EIG1R(:,:,1);
    EIG1R2 = EIG1R(:,:,2);
    EIG1R1(~Isel) = NaN;
    EIG1R2(~Isel) = NaN;
    EIG1R(:,:,1) = EIG1R1;
    EIG1R(:,:,2) = EIG1R2;
    % Eigenvalue 2, red:
    EIG2R1 = EIG2R(:,:,1);
    EIG2R2 = EIG2R(:,:,2);
    EIG2R1(~Isel) = NaN;
    EIG2R2(~Isel) = NaN;
    EIG2R(:,:,1) = EIG2R1;
    EIG2R(:,:,2) = EIG2R2;
    % Partial grid:
    X(~Isel) = NaN;
    Y(~Isel) = NaN;
end

hold on;

if Options.VisualizeGridHSA
    plot(X,Y,'g',X',Y','g','HandleVisibility','off');
end

quiver(X,Y,...
    visScale*EIG1B(:,:,1),visScale*EIG1B(:,:,2),0,strLineB,...
    'HandleVisibility','off','LineWidth',LW);   % (0: disable automatic scaling)
quiver(X,Y,...
    -visScale*EIG1B(:,:,1),-visScale*EIG1B(:,:,2),0,strLineB,...
    'HandleVisibility','off','LineWidth',LW);
quiver(X-visScale*EIG1R(:,:,1),Y-visScale*EIG1R(:,:,2),...
    visScale*EIG1R(:,:,1),visScale*EIG1R(:,:,2),...
    0,strLineR,'HandleVisibility','off','LineWidth',LW);
quiver(X+visScale*EIG1R(:,:,1),Y+visScale*EIG1R(:,:,2),...
    -visScale*EIG1R(:,:,1),-visScale*EIG1R(:,:,2),...
    0,strLineR,'HandleVisibility','off','LineWidth',LW);
quiver(X,Y,...
    visScale*EIG2B(:,:,1),visScale*EIG2B(:,:,2),...
    0,strLineB,'HandleVisibility','off','LineWidth',LW);
quiver(X,Y,...
    -visScale*EIG2B(:,:,1),-visScale*EIG2B(:,:,2),...
    0,strLineB,'HandleVisibility','off','LineWidth',LW);
quiver(X-visScale*EIG2R(:,:,1),Y-visScale*EIG2R(:,:,2),...
    visScale*EIG2R(:,:,1),visScale*EIG2R(:,:,2),...
    0,strLineR,'HandleVisibility','off','LineWidth',LW);
quiver(X+visScale*EIG2R(:,:,1),Y+visScale*EIG2R(:,:,2),...
    -visScale*EIG2R(:,:,1),-visScale*EIG2R(:,:,2),...
    0,strLineR,'HandleVisibility','off','LineWidth',LW);

function contourMod(X,Y,F,Options)
% Modified contour function 
if ~Options.SmoothContour
    [cc,hh] = contourf(X,Y,F);
    clabel(cc,hh);
else
    Fs = smooth2a(F,2,2);
    surf(X,Y,-99*ones(size(X)),Fs);
    view(2);
    shading interp;
end